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We investigate on the procedure of extracting a "spectral density" from mixed QM/MM calculations and 
employing it in open quantum systems models. In particular, we study the connection between the energy gap 
correlation function extracted from ground state QM/MM and the bath spectral density used as input in open 
quantum system approaches. We introduce a simple model which can give intuition on when the ground state 
QM/MM propagation will give the correct energy gap. We also discuss the role of higher order correlators of 
the energy-gap fluctuations which can provide useful information on the bath. Further, various semiclassical 
corrections to the spectral density, are applied and investigated. Finally, we apply our considerations to the 
photosynthetic Fenna-Matthews-Olson complex. For this system, our results suggest the use of the Harmonic 
prefactor for the spectral density rather than the Standard one, which was employed in the simulations of the 
system carried out to date. 



I. INTRODUCTION 

In the study of the dynamics of large systems such as pho- 
tosynthetic complexes, reduced models, which provide infor- 
mation on a small set of system degrees of freedom at the 
price of tracing out the rest of the bath degrees of freedom, 
have become very popular. Amongst these methods, which 
are open quantum systems approaches, one can find various 
quantum master equations [1-22] and Stochastic Schrodinger 
Equations [23, 24], which often rely on describing the system- 
bath interaction through a two-time "bath correlation func- 
tion" or a "bath spectral density". The reason for the quotes is 
that there is no unique way of obtaining these two quantities 
if one resorts to quasiclassical theories. In 1644, John Milton 
argued for the fact that "Reason is but choosing" [25], which 
inspires the philosophy behind this paper. Collecting informa- 
tion about the problem at hand may allow for better choices 
amongst the possible alternatives of the methods available to 
obtain these quantities. 

As a case study, in the present work, we consider the 
Fenna-Matthews-Olson (FMO) light-harvesting pigment pro- 
tein complex found in green sulfur bacteria. For this system, 
recent efforts have been undertaken to extract the bath spec- 
tral densities from mixed quantum-classical calculations [26- 
28]. The FMO complex has a trimeric structure, where each 
monomer contains, within the protein scaffolding, eight bac- 
teriochlorophyll (BChl) molecules, which can transport elec- 
tronic excitation energy. Up to recently, it was thought that 
only seven of the BChl's actually were present and most of the 
previous studies have focused on that case. Shim's results [27] 
which we employ in this work, are indeed based on the case 
of one monomer with seven BChl molecules. Experimentally, 
it has been possible to extract a spectral density which is an 
average over the BChls [29, 30]. However, one can expect that 
each BChl has a different spectral density, due to its specific 
protein environment. 
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A straightforward theoretical approach to obtain the spec- 
tral densities from a microscopic description is a mixed quan- 
tum mechanics/classical mechanics (QM/MM) model [31]. In 
this approach, the nuclear degrees of freedom are treated clas- 
sically and the relevant system quantities are calculated quan- 
tum mechanically. Then, from the microscopic description, 
spectral densities and correlation functions can be extracted 
and employed in the reduced models. 

A specific QM/MM approach, which has become popular 
in recent years in the context of photosynthetic complexes 
[26, 27, 32, 33] and has been employed for FMO [26, 27], con- 
sists in propagating the nuclei in the ground electronic state 
of the FMO complex, thus the change in the classical forces 
due to excitation of the BChls is ignored. The bath correla- 
tion function and spectral densities are then extracted from the 
energy gap trajectories, i.e. the electronic transition energies 
which depend on the time dependent nuclear configuration. 
This transition energy is calculated using quantum chemistry, 
for example TDFT [27] or semi-empirical approaches [26]. 
One thus obtains a time dependent energy gap two-time cor- 
relation function. 

In this work, we shed light on the connection between the 
mixed QM/MM gap correlation function and the open quan- 
tum system bath correlation function using a simple model. 
The mixed QM/MM gap correlation function is real. How- 
ever, in general, the full quantum correlation function will 
have an imaginary part. We employ different semiclassical a 
posteriori corrections to recover this part and compare the re- 
sulting spectral densities. Much work has been carried out on 
these a posteriori semiclassical corrections [34-39], but the 
question of which approximation is best remains open. To- 
wards answering this question, we show that a simple model 
of shifted harmonic Born-Oppenheimer surfaces leads to two 
of the semiclassical a posteriori corrections, each obtained 
with a different phase space probability distribution. We thus 
establish the link to a microscopic picture. This model of 
shifted harmonic potential surfaces is of particular interest, 
since the spectral density used in the open quantum system 
approaches emerges from such a description. 

Finally, we will investigate whether the results of the 
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QM/MM fulfill the requirements needed to employ the ex- 
tracted spectral density in open quantum system methods. 
These methods rely on the validity of assumptions such as lin- 
ear coupling of the system to the bath and a bath of harmonic 
oscillators. We attempt to investigate whether these assump- 
tions are valid in our case by evaluating higher order correla- 
tors of the energy gap time traces. In addition, we compare the 
spectral densities obtained at different temperatures. In most 
open quantum system approaches, when the harmonic bath 
approximation is employed, the spectral density is tempera- 
ture independent. Thus, one can use this invariance as a crite- 
ria for choosing which of the applied a posteriori corrections 
is most reasonable. In particular, our findings suggest that the 
best a posteriori semiclassical approximation for FMO is the 
Harmonic [40, 41] correction rather than the Standard [42^4] 
one, which has so far been employed in the context of the sim- 
ulation of exciton dynamics in photosynthetic complexes [26- 
28, 33]. Together, these aspects provide a clearer microscopic 
picture of the complex approximations involved in combining 
ground state QM/MM and open quantum system approaches. 

The paper is structured as follows: we begin by introducing 
the general quantum two-time correlation function in Section 
II. We introduce its time symmetries and its Fourier transform 
and subsequently we define the spectral density. A brief sum- 
mary of the general a posteriori semiclassical approximations 
to the quantum Fourier transform of the correlator from the 
classical Fourier transform is given in Section II C. In Sec- 
tion III, we introduce the concept of an energy gap correlation 
function for two-level systems as models for molecules cou- 
pled to a bath and show how this leads to a quantum bath 
correlation function and spectral density which are consistent 
with the open quantum system approach. In Section IV, we 
show that one can introduce a microscopic model which leads 
to some of the same prefactors described in the general case 
in Section II C. Finally, we investigate the conditions of lin- 
ear system-bath coupling and harmonic bath in Section V. 
In particular, we evaluate high-order multi-time correlation 
functions for the bath. These considerations are applied to 
our specific QM/MM calculations for FMO in Section VI. We 
conclude in Section VII by summarizing our findings. 



II. THE QUANTUM CORRELATION FUNCTION AND 
THE SPECTRAL DENSITY 

In this section, we introduce the definition of the quantum 
two-time bath correlation function. The generic Hamiltonian 
of a system coupled to a bath, in the absence of external fields, 
can be expressed as 

H = Hs (q,p) + i?B (Q,P) + ^SB (q,P, Q, P) , (D 

where iJg is the system Hamiltonian, H-q is the bath Hamil- 
tonian, iJsB is the system-bath Hamiltonian. In addition, 
(q,p) = {l]iP]) and (Q,P) = {Qk,Pk), indicate the gen- 
eralized multidimensional conjugated coordinates for the sys- 
tem and the bath respectively. The indexes j = 1, / and 
k = 1, F run over the system (/) and bath (F) degrees 



of freedom respectively. The system-bath Hamiltonian can be 
written as a function of the system. A, and bath, B, operators: 

HsB (q,P, Q, P) - ^ 1,„ (q,p) ® B,n (Q, P) • (2) 

m 

The influence of the bath on the system can be described 
by time-correlation functions. We will mostly focus on the 
two-time bath correlation function VII, which is defined as 

Cnm{t,t') = trB{B„(t,Q,P)A„(t',Q,P)/5B}. (3) 

Here, the time dependent bath operators are expressed as 
B™(t,Q,P) = e''^^^/^B„,{Cl,T?)e-'^^^/^. If the bath 
density matrix /5b is a stationary state of the reservoir, 
i.e. [-/Jb, Pb] — 0, then the correlator will be homogeneous in 
time [45], and depend only on the difference in time t = t = 
t', namely, 

C„™(t, O -C^nm(T) -trB{B„(r,Q,P)B„(0,Q,P)/5B}. 

(4) 

In the following we will be interested only in the n = m corre- 
lators, which we will indicate as C{t), dropping the subscript 
notation for simplicity. Finally, we define the equilibrium bath 
density matrix operator 



where 13 — 1/ {k-oT) and T is the temperature. In section V, 
we will briefly discuss higher order correlators. 

A. Fourier transform of the time correlation function and 
symmetries of tlie correlator 

The correlator defined in the previous section is in general 
complex and one can show, see e.g. [42, 46], that it has the 
following symmetries with respect to time, 

C{-t) = C* (t) = C{t - i/3h). (6) 

These symmetries result from Eqs 4 and 5. In addition, by 
looking at the expansion of the Fourier transform of C(<) with 
respect to h and using its second time symmetry in Eq. 6, one 
sees [35, 47] that the real and imaginary parts of C{t) are 
related by 

Im{C(<)} = tan(^^|)Re{C(i)}. (7) 

It is also useful to note that the real part of C(t) is sym- 
metric and the imaginary part is antisymmetric with respect 
to time. This connection is the starting point of many of the 
semiclassical approximations to the quantum correlation func- 
tion as described in Appendix A. 

We then define G(a;) the Fourier transform of the time cor- 
relation function 

/oo 
e"^'C{t)dt. (8) 
-OO 
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The function G{uj) is in general, temperature dependent, 
real and positive. In this work, we will refer to it as the 
Temperature-Dependent Coupling Density (TDCD). It will be 
convenient to split G{uj) into a symmetric and antisymmet- 
ric component which originate respectively from the real and 
imaginary parts of C(i), 

G(w) = Gsym(w) + Gasym(w), (9) 

with 

G,yM^\{G{Lo)+G{~uj)) (10) 

GasymM = ^(GH-G(-t^)). (11) 

In the definitions Eq. 9-11 we have followed the convention 
of Ref. [35]. Note that in the literature there exist other def- 
initions. E.g. the corresponding equations in Ref. [48], differ 
by a factor of 2 from the ones used here [49]. The detailed- 
balance condition, which follows directly from the second 
time symmetry in Eq. 6, implies that the overall TDCD is re- 
lated to its symmetric part by 

G{u:) = (12) 
= (1 - coth {phuj/2)) Gsy„(c^). (13) 
Alternatively, one can employ the asymmetric part to obtain, 

G{U) = ^_A_Gasy„.(^) (14) 
= (1 + coth [finw/2)) Ga.y,„(^). (15) 

We have expressed both equalities as a function of hyperbolic 
functions which will later be found in the various semiclassi- 
cal expressions of G(a;). 

B. The spectral density 

Another quantity which is often of interest is the so-called 
"spectral density". There are different definitions of spectral 
density in the literature. For later use in Section 11 C, we de- 
fine the spectral density j (w) as the positive frequency part of 

Gasym(^)[50] 

jH = Gasy„,H/ (27r) ; L^>0. (16) 

The spectral density j(a;) describes the coupling of the sys- 
tem to the bath for a given frequency. We extend the spectral 
density j(a;) to negative frequencies by defining 

J{u:) = Ga.y„.(a.) = I _27r . j (-co) ; c < " ^^^^ 

Using Eq. 14 and the definition of G{uj), Eq. 8, one can 
express the correlation function as a function of Gasym(w), 

1 f°° 

C{t) = — / dwe"*"* {coth{/3huj/2) + 1) J{u). (18) 

271- J-oo 



C. General semiclassical a posteriori approximations 

For systems of more than a few degrees of freedom, and in 
general, it is difficult to calculate the exact correlation func- 
tion, and therefore its Fourier transform, by using a fully quan- 
tum mechanical treatment. However, using classical mechan- 
ics one can obtain its classical counterpart with much less ef- 
fort. Therefore, it is common to attempt to construct the quan- 
tum spectral density from the classical one. 

Here, we define the fully-classical correlation function as 
the classical ^ ^ limit of Eq. 3, 

G'=\i) = J dQdPB (t,Q,P)B(0,Q,P)W(Q,P). 

(19) 

The classical bath phase-space density is defined as 

~fSHs(Q,P) 

Here, the quantum bath operators B in Eq. 3 have been sub- 
stituted by classical functions of the phase space variables 
B {t, Q, P). Also, the bath density matrix is replaced by the 
normalized equilibrium classical bath Boltzmann distribution. 
The classical TDCD is defined as 

/oo 
e*"*G'='(i)dt. (21) 
-oo 

Note that C'^^{t) is a real and symmetric function in con- 
trast to its quantum counterpart. This is also the case of the 
mixed QM/MM simulations employed for FMO in [26, 27], 
in which the resulting correlation function is real and none of 
the imaginary part is recovered. 

It is now desirable to be able re-construct, at least partially, 
the exact quantum spectral density from the classical one, 
through a simple description. Ideally, such a correction should 
be applied a posteriori and should not require extensive addi- 
tional computation. Much work has been carried out in this 
direction, see e.g. [34-39]. As described in Ref. [35], one can 
define various semiclassical approximations to the full quan- 
tum mechanical G(w) starting from its classical counterpart 
G'^'(w). We report each of these approximations in Table 1, 
second column, and a brief discussion can be found in Ap- 
pendix A. 

These corrections all originate from expansions in h and use 
of the symmetry properties of the two-time correlation func- 
tion and its Fourier transform. Note that if one expands the 
quantum correlator C{t) in powers of h, the first term is real 
and symmetric and corresponds to G'^^{t). The assumption 
that C{t) — G'^'(i), which leads to the standard approxima- 
tion, is in general not correct. In fact, since both of the corre- 
lation functions are obtained after thermal averaging, we see 
that they must differ at least by their respective partition func- 
tions. 

The specific functional form of each of these prefactors for 
G{uj), is plotted in Fig. 1 as a function of the temperature 
scaled frequency — Phu. We see that at low frequencies. 
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Table I. Column two: Various expressions for obtaining a semiclassical temperature-dependent coupling density TDCD G{(jj) from the 
classical G"^'(aj) as discussed in, e.g. [35]. Column three: Expressions for obtaining the semiclassical asymmetric TDCD J{uj) from the 
classical G'^'(ctj). These follow from the expressions in column two and from detailed balance (Eq. 14). 

Method Expression for G{lo) Expression for J (uj) — d^ymi^) 

Standard [42-44] G="(a;) = j:^^^ G"'(a;) J="'(a;) = tanh (^) G"'(a;) 

Harmonic [40-42] G''™(a;) = l_f-%nu, = ^G"'(a;) 

Schofield [51] G="'^°(a;) = ef'''"^^G''\uj) J="'^°(a;) = sinh(^)G"'(a;) 

Egelstaff [52] G'=s"'(a;) = e'^''"'/^ e*"*G"' (^f^ + (/3;i/2)2j dt r'^'=\uj) = sinh(^)J- [G"'(\/i' + {Pf^/'^Y)] {'^) 

Schofield-Harmonic [35] G=-''(a;) = e^^'^/^'y^^^f^^^G^H'^) = y'^sinh(^)G"'(a;) 




Figure 1 . Prefactors for the temperature-dependent coupling density 
(TDCD) G{u}) (Tab. I, column two) plotted as a function of tUb = 
ujhp. One sees that for huj/kBT < 1 they coincide. The prefactors 
for the asymmetric J (to) = G^°^™(a;) (Tab. I, column three) follow 
the same trend with the difference that J(0) — 0. 



cjb < 1 (i.e. hu> < ksT) all approximations give nearly iden- 
tical results and give the same value for = 0. 

The various approximations to the spectral density can 
straightforwardly be derived from those of G(w) by using 
Eq. 14. The resulting expressions are reported in column three 
of Table I and the prefactors follow the same trend as those for 
G(ti;) as a function of frequency. 

Now, given all the functional forms described above, the 
question is how to reason and choose the most appropriate 
one. For the FMO complex, it is unclear at first sight which 
one would be the best. In Section IV, we will investigate a 
model to elucidate the origin of these prefactors. This will 
help to discriminate between these corrections. In Section VI, 
we will apply all of the corrections listed in Table I to our 
energy gap traces and discuss the differences between each 
approach. 



III. ENERGY GAP CORRELATION FUNCTION FOR A 
SIMPLE MODEL 

In the mixed QM/MM calculations for photosynthetic sys- 
tems [26, 27, 33], the nuclear trajectories are propagated in the 
electronic ground state using MD with short time steps. For 
a set of longer times steps within these trajectories, the elec- 
tronic transition energies of the BChl molecules are computed 
using an electronic structure calculation method. Because it 
is computationally costly to calculate the electronic states for 
the full set of seven/eight coupled BChls simultaneously [26], 
the system was divided into seven/eight subsystems for which 
the electronic states were calculated separately. Thus, in these 
calculations no excited state interactions are included explic- 
itly [53]. The Hamiltonian of the coupled BChls is then writ- 
ten as H = Y.n=i Hn + Y.n<ni ^here i?„ denotes the 
Hamiltonian of BChl n and Kimis the Coloumb (transition 
dipole-dipole) interaction between them. To establish a con- 
nection to the open quantum system approach, each BChl is 
treated as an electronic two level system. These two-level sys- 
tems and the electronic interaction between them are taken 
to be the system part. The coupling to internal nuclear de- 
grees of freedom an the surrounding protein will then lead to 
fluctuations of these quantities in time (for more details see 
e.g. [33]). From the time dependence of the transition energy 
between electronic ground and excited state for each BChl, a 
classical ground-excited state energy-gap correlation function 
can be obtained. In turn, spectral densities can be extracted 
from the energy-gap correlation functions. 

The gap correlation function, as obtained from the IVID sim- 
ulations, is a quantity which up to the previous section, has 
not been connected to the open quantum system approach de- 
scribed in Sec. II. To this end, in this Section, we will explore 
a simple model with Born-Oppenheimer (BO) surfaces which 
can clarify the connection. 



Quantum correlation function and energy gap correlation 
function for a molecule 



Lets us begin by considering a single molecule (BChl) 
treated in the Born-Oppenheimer approximation. The 
molecule is modeled as a two-level system with an electronic 
adiabatic ground \g) and excited |e) state. This situation is 
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Figure 2. Sketch of two generic Bom-Oppenheimer (BO) surfaces. 
This sketch indicates the model which we consider in Sec. Ill A to 
define the energy gap correlation function in the context of an arbi- 
trary molecule. 

depicted in Fig. 2. We can think of the BO-surfaces as having 
the dependence of the environment (protein and other BChls) 
akeady included, ignoring however the resonant dipole-dipole 
interaction. The approximation of two levels is reasonable in 
the limit where the next excited state is very far in energy 
space from the first. Usually, non-adiabatic couplings can 
be also neglected, as chosen in our calculations. Given this 
model, we investigate how the general correlation function, 
Eq. 3, is related to the energy gap correlation function. 
We write the full Hamiltonian formally as 

H = ffg(Q,P) \g){g\ + 7?,(Q,P) |e)(e| , (22) 

where Hg{Q) and HdQ) are the nuclear Hamiltonians for 
the ground and excited state in the BO approximation. In mass 

scaled coordinates (^Qj = ^nijqj] Pj — the Hamilto- 

F p2 

mans can be expressed as Hg{Q,P) + Vsr{Q) 

and HdQ,P) = Hg{Q,P) + Aog(Q) where Vg(Q) denotes 
the grounds state potential energy surface. For later purpose, 
we have expressed the excited state nuclear Hamiltonian with 
respect to the ground state potential by introducing the energy 
gap operator, 

Aeg(Q) = HoiQ.P) - H^iQ, P) (23) 

= ;iCJeg + Ao + Fc(Q)-V^g(Q). 

This operator quantifies the energy difference between the ex- 
cited state and the ground state surface. A coordinate indepen- 
dent constant energy difference huteg + Aq has been explicitly 
written down, so that the remaining partly (Q) — V'g(Q) does 
not contain any coordinate independent contributions. This 
notation and the division into two parts hui^g and Aq will be- 
come clear later 



The total Hamiltonian can be rewritten as 

H = K,-i + (fic^eg + Ao) |e)(e| + A |e)(e| , (24) 
where we have defined the reduced gap operator A = Acg — 

/iWog — Aq. 

To establish a connection to the open quantum system 
model, as presented in Sec. II, we choose 

- 4(Q.P) (25) 
i?SB - A(Q) |e)(e| (26) 
= (fi^cg + Ao) |e)(e| , (27) 

where we have set the energy of the electronic ground state 
\g) to zero. From the form of H^b we identify the system 
operator Ac = |e)(e| and the bath operator B — Aog(Q). We 
can now define the usual bath correlation function as 

C{t) = trB{B(t)B(0)pB} - trB {A(i)A(0)pB} , (28) 

where we have dropped the dependence on bath coordinates in 
the notation for simplicity. A can be thought of as a "gap" op- 
erator, that is, as a measure of the energy difference between 
the ground and excited state at a given nuclear configuration. 
From now on we will indicate reduced gap correlation func- 
tions as 

a(0 =trB{A(t)A(0)/5B}, (29) 

to distinguish them from the general bath correlation function 
C{t). Eq. 29 corresponds to the full quantum gap correlation 
function that one would obtain, e.g., from a quantum simu- 
lation on the FMO complex, considering only two electronic 
levels per molecule and after including the protein environ- 
ment. 



B. Quantum correlation function and energy gap correlation 
function for harmonic surfaces 

In most of the open quantum system approaches used to 
describe the FMO complex, the bath is taken as an (infi- 
nite) set of harmonic oscillators for the environment of each 
BChl. Each oscillator coordinate is then assumed to be lin- 
early coupled to the electronic excitation of the BChls, i.e. 
^^SB = |e)(e| ® i^jQi where Kj is a coupling constant 
defined in Tab. II. To establish the connection between the 
reduced gap operator and this system-bath interaction, we 
now consider identical shifted harmonic potential surfaces, as 
sketched in Fig. 3. The quantities defined in the general case 
in the previous Section III A become. 



Fe(Q,P) = hw,^ + \ Y^iPf + ^]{Qj - 5Qjf) 

j 

= 7?g(Q,P) + Aeg(Q), (30) 
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Figure 3. Shifted liarmonic identical Born-Oppenheimer surfaces. 
This model is the one employed in Sec. IIIB to derive classical and 
semiclassical expressions for the Fourier transform of the bath corre- 
lation function G{uj) and for the spectral density. 



with 



Aeg(Q) = fit^eg + 2^PQ'' - ^^(^^'^g.OQ,, (31) 

i 3 



in Eq. 17 with 

F 

jH (34) 

In consistency with our definition of the Fourier transform, 
Eqs. 8 and 18, a factor of 2n has been introduced and appears 
in Eqs. 33 and 34. Note that a{t) corresponds to the C{t) of 
Eq. 3. 

To establish a connection to the classical correlator, which 
is real and symmetric, we note that can be obtained from 
the real part of a{t) via 

= ^tanh^^^^ ^ Rc{a{t)}cos.{ut)dt. (35) 

This is the expression (up to the constant prefactor 2/h ) used 
in Refs. [27, 33], to obtain spectral densities. When using this 
expression to obtain the spectral density from QM/MM simu- 
lations one often assumes that C^'(t) « Re{a(t)}, following 
the Standard approximation (see. Appendix A). Then, after a 
Fourier transform and use of symmetry relations for G{ll:) one 
finds the following expression, 

j{uj) = i- tanh(^^) G^\uj) ; O 0. (36) 

Note that in this section, Eqs 33, 35 and 36, differ by a 
factor of 1/2 from those in Ref. [48] due to our definition in 
Eq. 11. 



where ilj is the frequency of the j-th oscillator Here, 
the constant shift Ao introduced in Eq. 23 corresponds to 
the frequently empoyed reorganization energy Aq = A^ = 
\^]5Q^- Using Qj = ^h/{2nj){a] + a^) and Pj = 
iy/h^lj /2{a^j — Gj), Eqs. 30 can also be written as 



He. — He: 



Here Oj and aj are the bosonic creation and annihilation op- 
erators coiTesponding to the ground state Hamiltonian. We 
have also introduced the so-called Huang-Rhys factor [54]: 
Xj = ^SQ^. Note that the total Hamiltonian is now in the 
standard form of an open quantum system model, as in Eq. 1, 
with the relevant quantities given in Table II. 

This model for a few vibrational modes of the chro- 
mophores, has been successfully employed to describe the op- 
tical properties of molecular aggregates [55-58]. 

From the expression of the energy gap operator shown in 
Table II, one obtains the quantum two-time bath correlation 
function [48] 



hide 




a(t) = 2 / J[io) cotn cos [ujt) — isin(wi) 

2vr Jo L V 2 / 

(33) 

where J(a;) = 27r {<d{uj)j{uj) — j(— a;)8(— oj)) , was defined 



IV. CLASSICAL AND SEMICLASSICAL LIMITS OF THE 
CORRELATORS AND SPECTRAL DENSITIES FOR 
HARMONIC SURFACES 

As outlined in the previous section, the harmonic model al- 
lows for a simple analytic solution in the quantum mechanical 
case. Now we will show that the system also has a solution in 
the classical case. In particular, in this section, we will intro- 
duce a model to construct exact relations between the classi- 
cal gap-correlation and the quantum one. To this end, we will 
consider classical dynamics in the ground state BO potentials 
within an initial value representation of the initial state which 
is consistent with the mixed QM/MM approach. For each 
initial value, we calculate a trajectory and the corresponding 
reduced classical energy gap between the two surfaces, i.e. 
A(Q(i), P(^)). We then average over many trajectories. 



A. Classical equations of motion 

The classical equation of motion of the j-th harmonic 
bath coordinate is Qj + il^Qj — 0. Solving this dif- 
ferential equation with the initial condition (Qjo^Pjo) — 
{Qj{t = 0), Pj{t — 0)) yields the time dependent coordinate 
trajectories 



]j{t; Qjo, Pjo) 
]jo cos(iljt) + 



PjO 



sin(f7j<). 



(37) 
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Table II. Expressions of the system bath quantities for the case of two Bom-Oppenheimer harmonic surfaces as sketched in Fig. 3. 



Quantity 



Expression 



System Hamiltonian 
System-bath Hamiltonian 
Bath Hamiltonian 

Energy gap operator 



Hs = {hojcg + \r) \e){e 
HsB = \e){e\ A{Q) 



Reduced energy gap operator A(Q) — Acg(Q) — {huJag + Ar) 
Reorganization energy Ao = Ar = J]] . ^Q,^SQ^ 



Coupling constant 
Huang-Rhys factor 
Unitless constant 



For each trajectory, the energy gap is then expressed as 

A(i) = A(t;Qjo,Pjo) 



with yV]'°'*^(Qjo, Pjo) = ^e^^^-^i'o+^^^Qjo)^ and it is nor- 
malized to one, i.e.,/ dPjodQjoyV^°^^^{Qjo, P]o) = 1- Note 
that (rj2^Qj)2 = 2hXjn^y Using Eq. 39 and the Boltzmann 



= — ^(r2|(5(5j) ^Qjo cos{VLjt) + sin(fijt) j distribution for initial positions and momenta, we obtain 



(38) 



where the parametric dependence of Qj and A on the initial 
conditions (QjOi ^jo) has been explicitly indicated. 



B. Energy gap correlator 

The evaluation of the reduced gap correlation function, 
Eq. 29, in the classical limit, results in the following expres- 



= E J dPodQoWiQo^Po) 

jk 



Ait;Q,o,Pjo)A{0;Qko,Pjo) 



(39) 



where yV(Qo, Po) is the initial distribution and dPorfQo de- 
notes the set of all coordinates, i.e. dQo = dQiQ ■ ■ ■ dQuo- 
For harmonic potential surfaces, Eq 19, is time-evolved fol- 
lowing Eq. (38). In this Section, we will investigate two dif- 
ferent choices for the initial distribution, namely a Boltzmann 
distribution, as in Ref. 21, and a Wigner distribution which 
resembles the quantum thermal state. We will refer to the two 
cases as the classical limit and the semi-classical limit, respec- 
tively. 



C. Classical and semiclassical correlation functions 

1. Classical limit 

To obtain the classical limit of the correlator, we choose the 
Boltzmann distribution for the initial coordinates which cor- 
responds to a purely classical thermal state. The distribution 
is defined as follows 



boltz 



(Qo,Po)=n Wi°"'(QjO'^:'0 



(40) 







(41) 



Here we have introduced the definition Qj = Mlj/ [kBT) . 
This quantitiy is also reported in Tab. II. 



2. Semiclassical limit 

In order to obtain the semiclassical limit, we take the quan- 
tum Wigner distribution for the initial coordinates and use it 
in Eq. 39. The Wigner distribution is given by 

W-'s(Qg^ p„) = -Q wfs(Q,o, P,o), (42) 



where we have used the compact notation 



W7'S(Q,o,-P,o) = 2tanh| 







(43) 



The normalization of the Wigner distribution is chosen such 
that S ^^^Wf'^{QjQ,P.jo) = 1. The resulting expres- 
sion of the energy gap correlation function is 

awig(i) = > ]{hVtjYX.j cos{fljt) coth( ^ ) . (44) 



E 



D. Classical and semiclassical spectral densities 

The Fourier transform of the correlator a{t), Eq. 8 yields 
the TDCD G(w). After applying it to Eq. 41 and to 44 we 
obtain 



(45) 
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and for the Wigner distribution 



Gwig (uj) — TT ^ coth 



hid 



(46) 



When only terms up to first order in Q are significant, as in 
the case of the Harmonic surfaces in the hnear system bath 
coupHng Hmit, Tab. II, we can write the two-time correlation 
function as 



HereK, =;ia;yX"asinTab. II. Now,usingEq. 34wehave 0) / A(i) A(0) ■) ^ {Q{t)Q]{^))- (51) 



Gboitz(w) = 27r ^^^"^ j(a;) 



Gwig(w) 



fujj 
2-K cothl 



(47) 
(48) 



One can then obtain the corresponding spectral densities 

We see that the semiclassical Wigner distribution yields the 
same prefactor as for the Standard approximation described in 
Section II C while the Boltzmann distribution gives the same 
prefactor as the Harmonic approximation, also described in 
Section II C. 



V. MODELS FOR SYSTEM-BATH COUPLING - HIGHER 
ORDER CORRELATORS 

As discussed in the introduction, there has been a lot of 
interest in modeling the exciton dynamics of the FMO com- 
plex using open quantum system approaches. These usually 
require as input a bath two-time correlation function or (equiv- 
alently) a spectral density and they rely on the assumption of 
linear coupling to the bath and on a bath described by har- 
monic oscillators [59]. 

In the previous Section IIIB, we have discussed that this 
model corresponds to shifted adiabatic BO surfaces of iden- 
tical curvature. We have shown that in this case, the energy 
gap two-time correlation function for a classical ground-state 
propagation is directly proportional to the quantum one and 
we have extracted the appropriate (frequency dependent) pro- 
portionality constant. For other shapes of the potential sur- 
faces involved, one will in general obtain different propor- 
tionality constants, although the delta-peaks of the spectral 
densities can be located at the same energies (the positions 
are determined by the shape of the ground state potential) . 

It is not clear, a priori, if the approximation of shifted 
harmonic surfaces (or equivalently linear coupling to a har- 
monic bath) is a good one for the system under consideration. 
To gain some insight on this question, from an analysis of 
QM/MM trajectories, one possibility is to consider higher or- 
der correlators. If the approximation of linearly coupled har- 
monic oscillators is inadequate, one expects that higher order 
correlators will have a significant relative weight. 

We proceed to discuss some properties of correlations of the 
bath gap operator, Eq. 23. The energy gap operators can be 
described by a function of the bath coordinates and expanded 
in terms of these as 



(50) 



Here, we have excluded the zeroth-order term which corre- 
sponds, e.g. to a reorganization energy, and is usually renor- 
malized into the system Hamiltonian. The angular brackets 
(...) = tre {..., pb} indicate thermal averaging over the bath 
degrees of freedom. Similarly, the three-time correlation func- 
tion becomes 



(i',i,0) = (A(i')A(i)A(0) 



= ^d'^ef ^ef ^ mmAmm) ■ (52) 

In the case of a harmonic bath, the three-time correlation func- 
tion will vanish, and in general any odd permutation of the 
harmonic bath coordinates will vanish. 

However, if one considers the case where one retains the 
second order term in Eq. 50, the two-time correlator will be- 
come: 

a(i,0)= SyS;,,(Q,j.(i)QH(0)), (53) 

ijfci=0 

where we have defined S^j and Q^j (t) as 

'0 ;i-j = 

; J - A * ^ 
if ; i = A J ^ 




; « - J - 

; J = A z ^ 
; z = A 



Q,{t)-Q,{t) ;i,j^O 



Analogously the three-time correlator becomes 



ijklmn—O 



^ij^kl^rnn (^^^(^ ) Qfc/ (^) Qmn (0)) ■ 



(54) 

If the bath is harmonic, it is straightforward to show that 
all terms with an odd number of coordinate operators in the 
averages will vanish. Yet, we see that in general, unless the 
coupling to the bath coordinates is linear and the bath con- 
sists of Harmonic oscillators, the three-point correlator will 
not vanish. It may therefore be necessary to go beyond the 
simple description using only the two-time correlator. 



VL APPLICATION TO THE FMO COMPLEX 

In this section, we apply the approximations discussed in 
Section II C, to the energy gap trajectories obtained from the 
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mixed QM/MM simulations for the FMO complex of Prosthe- 
cochloris aestuarii as carried out recently by us in Ref. [27]. 
The nuclear trajectories were obtained by classical MD using 
the AMBER 99 force field. An isothermal-isobaric (NPT) en- 
semble was employed in the MD simulations. For the calcu- 
lation of the energy gap, snapshots of the nuclear coordinates 
were taken at every 4 fs. For each ground state configuration, 
the gap was obtained by computing the energy correspond- 
ing to the Qy transition of the BChl's using time-dependent 
time-dependent density functional theory with BLYP func- 
tional within the Tamm-Dancoff approximation. More details 
on the computation can be found in Ref. [27]. 



A. TDCD and spectral density from mixed QM/MM with a 
posteriori semiclassical corrections 

Using the energy gap trajectories obtained in Ref. [27], we 
evaluated the different semiclassical approximations as re- 
ported in Tab. 1. We denote the time-points at which the en- 
ergy gap is calculated by ti and the corresponding energy gap 
by Xi where i = Q . . . N — \ runs over the N the time-points. 
As in Ref. [27] we evaluated the correlator by using a discrete 
representation, which implements the fc-th element of the two- 
time correlator as 

N-k 

^^- (Ar-fc)^(2) E(-^-"-^) (^^+^-^) (55) 

where irS^^ — X]f ~ X)^/N is the variance and X is 
the mean. Here, one assumes that the N — k values Xi give 
a faithful initial distribution which reproduces the Boltzmann 
distribution. To minimize spurious effects in the Fourier trans- 
form, we multiplied the time trace by a Gaussian of variance 

^laussian = 0.09 • i^,, = 2.304 • 10^ fs^ with t„,ax = 1600 fs, 

the length of the correlation function (as reported in [27]). The 
Gaussian is normalized to have unitary area in frequency do- 
main [60] following our definition of the Fourier transform in 
Eq. 8. Next, we computed the different semiclassical quanti- 
ties of Table I using our initial time trace. 

In Figure 4, we show the temperature-dependent coupling 
densities TDCDs (as defined in Eq. 8), for site 1 of the FMO 
complex (site 1 at 77K and 300K) evaluated using the different 
approximations listed in Table 1 column two. We notice how, 
as expected, there are little differences between the approx- 
imations at low frequencies. Only at higher frequencies the 
TDCD differs significantly for each approach. The Egelstaff 
approximation incorrectly predicts a negative spectral density 
for low frequencies in this case and was therefore not shown 
in the plots. 

Similarly, in Fig. 5, we show the asymmetric component of 
the TDCD ( J(w) = G'asyin(t^)) for site 1 of the FMO complex 
at 77K and 300K evaluated with the different approximation 
methods. The trend is similar to what we saw for the TDCD. 
From the general definition of each semiclassical correction, 
it isn't clear which one is most accurate. To better reason on 
which one to choose we evaluated the three point correlator 
(Sec. VIC) to see if the approximation of a Harmonic bath is 



valid. The motivation for this procedure was given in Section 
V. 



B. Analysis of prefactors in terms of temperature dependence 
of the spectral density 

From our discussion in Section I, we recall that many open 
quantum system approaches rely on the assumption of linear 
coupling to a bath of harmonic oscillators. This leads to a 
temperature-independent spectral density. Inspection of Fig 
5 shows that the asymmetric TDCD (which is directly pro- 
portional to the spectral density) obtained from the QM/MM 
is not similar at different temperatures. This is more appar- 
ent at higher frequencies. The source of this temperature de- 
pendence could be due to the fact that the MD bath is not 
fully harmonic or that the coupling is not completely linear 
(it might also be possible that the MD is not fully converged). 
The reason for this temperature dependence is still not imme- 
diately clear to us, and we leave this for further studies. We 
would need to run longer QM/MM trajectories to improve the 
statistics and further check convergence of the distributions. 
Nevertheless, based on the most common assumption in the 
open quantum system approaches, we suggest that the relative 
temperature-independence is a criterium to select the best ap- 
proximation for the spectral density. In Fig. 6, we compare 
the asymmetric TDCD ( J(a;) = 27rj(aj) ; a; > 0) obtained us- 
ing the Standard and the Harmonic approximations. It is clear 
by comparing panel a) and b), that only the Harmonic cor- 
rection leads to a temperature-independent spectral density. 
In the bottom panels, c) and d) we compare the asymmetric 
TDCD averaged over the seven sites, with the experimental 
asymmetric TDCD. The asymmetric TDCD is obtained from 
the spectral density as defined in Eq. 17. Further, since in 
Refs [29, 30, 48], Eq. 11 is defined without the factor 1/2, 
the experimental spectral density is rescaled by 1 /2 for com- 
parison. From panel e) and f) of Fig. 6, we see that the Har- 
monic curve is in good agreement with experiment. Results 
for all sites at both temperatures are reported in the Supple- 
mentary information. Appendix B. The standard asymmetric 
component J(aj) = G'asym(i^) shown in Fig. 6 panel a), corre- 
sponds to the spectral density in Ref. [27] if one multiplies the 
Shim [61] resuh by 27r-tanh(a;/3/i/2)/tanh(j//3/i/2) to obtain 
J{ijj) (here v = a;/(27r)). The analysis of the temperature de- 
pendence of the prefactors for the spectral density suggests 
that the Harmonic approximation is the correct choice for the 
FMO complex rather than the Standard one, even though in 
the later case the agreement with the experimental FLN nar- 
rowing results is better. This analysis, is bolstered by the ev- 
idence provided by the model in Section IIIB, and the higher 
order correlator results, which will be presented in the next 
Section. 



C. Higher-order correlation function 

From the theory of discrete processes, similarly to Eq. 55, 
we see that the (fc, j)-th element of the three-time correlator 
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Figure 4. Temperature-dependent coupling densities G{uj) obtained with eacli of tlie Standard, Harmonic, Scliofield and Scliofield-Harmonic 
corrections as listed in Tab. I, column two. In particular in panel a) results are at 77K and in panel c) 300K. Panels b) and d) display the low 
frequency region of figures a) and c) respectively. These results are obtained after multiplication of the time series with a Gaussian window 
function as described in the Section VI A. 



is 

(iV-fc-j)m(3)'' 

{Xi — X) (Xi+k — X) {Xi-^-k+j — X) 

i=l 

(56) 

where m^'^) = {Xi — X)^ /N is the third moment and 
X is the mean. We compare the two correlators by dividing 
them by increasing powers of the standard deviation VmJ^. 
The results for site 1 of the FMO complex at 77 and 300K are 
reported in Fig. 7. For the two-time correlator. Fig. 7 panels 
a) and b), we see correlations up to at least 1000 time steps, 
while for the three-time correlator, panels c), d), e) and f), we 
see a rather flat profile with values about two orders of mag- 
nitude smaller than the largest value of the two-time correla- 
tions. This is observed for all sites and temperatures (Results 
for all sites can be found in the Supplementary information. 
Appendix B). 



This means that since we find a small and roughly flat three- 
time correlator, the linear coupling to a harmonic bath as- 
sumption is probably good. In fact, as described in Sec. V 
this case corresponds to linear coupling to the bath and Gaus- 
sian correlated bath operators. Of course, there may be other 
fortuitous cases in which the three-time correlator is roughly 
zero and the bath is not harmonic. Further, this comparison 
is based on the order of magnitude of the correlations, the 
three-time correlator is only much smaller It may be that 
for some modes of the system, certain frequencies, present in 
the three-time correlator's two dimensional Fourier transform 
give a more important contribution to the dynamics than other 
frequencies present in the spectral density. Nonetheless, the 
above result encourages the idea that the assumption of linear 
coupling and harmonic bath is valid. This, in turn, implies that 
one should use the Harmonic semiclassical correction in Sec. 
11 C, which is also consistent with the prefactor found in 111. 

One a final note, to confirm with certainty that the bath is 
Harmonic, one should evaluate higher order correlators, be- 
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Figure 5. Asymmetric component of temperature-dependent coupling density J (to) = G'asym(<^) ; obtained with each of the Standard, 
Harmonic, Schofield and Schofield-Harmonic corrections as listed in Tab. 1, column three. In particular in panel a) results are at 77K and in 
panel c) 300K. Panels b) and d) display the low frequency region of figures a) and c) respectively. These results are obtained after multiplication 
of the time series with a Gaussian window function as described in the Section VIA. Note that the spectral density, as defined in Eq. 17 can 
be obtained by dividing J{lj) by 2n. 



yond the three-time correlator. However, to obtain a statisti- 
cally relevant estimate, much longer time dependent energy 
gap trajectories, which are expensive in terms of the QM/MM 
propagation, would be required. Work in this direction is be- 
ing carried out in our groups. 



VII. CONCLUSIONS 

In this work, we have investigated the connection be- 
tween the gap correlation function extracted from ground state 
QM/MM and the bath spectral density used as input in many 
open quantum system approaches. 

One important point is that the classical bath correlation 
function is real while the quantum mechanical one is gener- 
ally complex. There exist several semiclassical a posteriori 
corrections which aim to fix this and we have employed them 
on our time traces to recover a part of the imaginary compo- 
nent. 



The discussed prefactors originate from general expansions 
in orders of h and do not include information on the specific 
type of system-bath coupling, etc. We have investigated two 
simple models and found that the prefactors obtained corre- 
spond to two of the general semiclassical expressions. Thus, 
we have linked the semiclassical limits with a microscopic po- 
tential energy surface picture. 

We have shown that the gap-correlation function extracted 
from ground state QM/MM only corresponds to the fully 
quantum excited state calculations in the case of shifted 
parabolas. This model for a few vibrational modes of the 
chromophores has been successfully used to describe the opti- 
cal properties of molecular aggregates. Including only a finite 
number of internal vibrations is probably a good approxima- 
tion for molecules in the gas phase or suprafluid Helium nan- 
odroplets [58]. However, for molecules in solution or when 
a protein environment is present it is no longer a good ap- 
proximation to include only only a few (undamped) modes. 
In particular, one has to take into account the interaction of 
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Figure 6. Panel a) comparison of tlie asymmetric component of temperature-dependent coupling density J{uj) = G'asym(i^) ; for site 1 for 
the Fenna-Matthews-Olson complex, obtained with the Standard approximation (Table I first line third column) at 77K and at 300K. Panel 
b) comparison of J(aj) obtained for site 1 with the Harmonic approximation (Table I second line third column) at 77K and at 300K. We 
see clearly that the Harmonic prefactor gives a roughly temperature independent J(w), while large differences are seen using the Standard 
prefactor. This consideration can also be made for plots in panels c) and d) where one sees J{ijj) averaged over all sites. In particular, Panel c) 
shows J{Ld) averaged over all seven sites calculated with the Standard approximation at 77 and 300K and the green curve corresponds to the 
experimental spectral density reseated by n (a factor 2-k comes from Eq. 14 and a factor 1/2 comes from Eq. 1 1 ) to obtain J{ijj) as defined 
in Eq. 17 [29, 30], Panel d) shows J{ijj) averaged over all seven sites calculated with the Harmonic approximation at 77 and 300K and again 
the green curve corresponds to the experimental spectral density [29, 30]. The agreement with the experimental (green) spectral density is 
better with the Harmonic approximation than with the Standard approximation. Panels e) and f) correspond to the same quantities as those of 
panels c) and d) in the low frequency region, here we note that both approximations are roug hly equivalent for ^ < 1 (e.g at T = 77 K for 

oj < 55 cm~^and at T = 300 K for a; < 200 cm~^). Note that the spectral density, as defined in Eq. 17 can be obtained by dividing J(aj) by 
27r. 
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Figure 7. Panel a): Two-time correlation function of the energy gap fluctuations of site 1 of the FMO complex, normalized by m'^' at 77K as 
in Eq. 55. Panel b): Two-time correlation function for site 1 at 300K. Panel c) Three-time correlation function of the energy gap fluctuations 

of site 1 of the FMO complex, as defined in Eq. 56 multiplied by m'^'/ mt^)^ at 77K. Panel d) Three-time correlation function for site 1 
at 300K. 
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the vibrations with the environment in addition to the direct 
interaction of the electronic excitation with the environment. 
For this general situation, it is no longer clear if the model of 
shifted harmonic potential surfaces is indeed a good descrip- 
tion of the system. 

Therefore, we have investigated whether the approxima- 
tion of harmonic bath and linear coupling is accurate for our 
QM/MM calculations for the FMO photosynthetic complex 
by computing the next higher order correlator beyond the two- 
time correlator. The three-point correlator seems to give a 
small contribution which, while not being conclusive, sug- 
gests to us that the Harmonic/linear coupling model is a good 
approximation. The analysis of the temperature dependence 
of prefactors for the spectral density also suggests that the 
Harmonic approximation is preferred to use for the FMO com- 
plex, and perhaps other photosynthetic complexes, rather than 
the Standard one. Having made these choices, the theoretical 
results are in reasonably good agreement with the experimen- 
tal spectral density. These result in a much better agreement 
than in our previous work, which underestimated the magni- 
tude of the spectral density [27] and than other QM/MM cal- 
culations [26] which overestimate the coupling to the bath by 
one order of magnitude. 

Finally, we have explained the link between bath correla- 
tion function and gap correlation function and found mod- 
els under which the gap correlation function can actually be 
viewed as a general open quantum system bath correlation 
function. 



ACKNOWLEDGMENTS 

We thank Dr Semion Saikin, for his useful advice and 
Gerhard Ritschel for going over the manuscript. A. E. ac- 
knowledges financial support from the DFG under Contract 
No. Hi 872/1-1. S. V. and A. A. G. acknowledge support 
from the Center for Excitonics, an Energy Frontier Research 
Center funded by the U.S. Department of Energy, Office of 
Science and Office of Basic Energy Sciences under Award 
Number DE-SC0001088 as well as support from the Defense 
Advanced Research Projects Agency under award number 
N66001-10-1-4060. S. V. and A. A.-G. acknowledge support 
from the Harvard Quantum Optics center. 



Appendix A: Semiclassical approximations to the quantum 
TDCD 

The first and most simple approximation is known as the 
standard approximation [42, 44, 62], and consists in assuming 
that C{t) ~ C"^'(i), i.e one assumes that the quantum correla- 
tion function is roughly equal to the zeroth order term in its h 
expansion. Then, taking the Fourier transform and using the 
appropriate Eq. 12 one finds. 



We see that this becomes a good approximation for small h 
or also, for small lu. 

The next approximation is known as the "harmonic" ap- 
proximation [40-42], which corresponds to taking the first or- 
der term in the expansion of Eq.7 with respect to h, and then 
using C{t) ~ C^'(t) . After Fourier Transformation and sub- 
stitution using Eq. 14, one arrives to the final expression Note 
that this approximation doesn't assume a harmonic bath, but 
turns out to be exact in the case of a system coupled to a bath 
of harmonic oscillators as shown in Ref. [35]. 

The third approximation is due to Schofield [51] and comes 
from a combination of the standard and harmonic approxima- 
tions, expanding to first order separately the imaginary and 
real parts of C{t) and Fourier transforming to obtain 

G''''^°{uj) = e'^'"^/^G''\uj). (A2) 

The fourth approximation is the Egelstaff approxima- 
tion [52] which originates from assuming that C{t) — 

C"^' ^ ■\/ {t(t + ij3h))^ and taking the Fourier transform to ar- 
rive at the final expression 

G'=s'='(w) = e'"'"/^ j e^^t^ci (^^2 + {ph/2)A dt. 

J — OO 

(A3) 

The same ansatz on C{t) can be used to obtain the Schofield 
approximation, by expanding to first order and subsequently 
Fourier transforming. 

Finally, the last approximation we wish to recall can be ob- 
tained by taking the logarithmic average of the Schofield and 
Harmonic approximations thus obtaining 

G^-(.)=e^-/^y^^lG-(.). (A4) 

From these derivations, it is clear that in general it is diffi- 
cult to identify the best approximation. The Egelstaff approx- 
imation doesn't always satisfy the fact that G{ijj) > while 
the other approximations do. Then, all apart from the Egel- 
staff have the constraint that G{0) = G"^'(0). As described 
in Ref. [40], one can also check that the Harmonic approx- 
imation is exact in the case of linear coupling to a bath of 
harmonic oscillators while the Egelstaff approximation is the 
most accurate for exponential coupling to harmonic coordi- 
nates. For the other expressions it is unclear which will be 
best for a generic bath. 

Appendix B: Supplementary information 

1. Two and tliree-time correlation functions for aU site of the 
fmo complex 

In this Section we report the two and three-time correlation 
functions for the energy gap fluctuations of the FMO com- 
plex as discussed in the text for all seven sites of the com- 
plex and for both temperatures, 77 and 300K. The functions 
are rescaled by increasing powers of the standard deviation 
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V m(2) as a means of comparison. We notice, by compar- 
ing panels a) and b) to c) and d) in Figures 8-14, that the 
three-time correlation function is always much smaller in am- 
plitude than the corresponding two-time autocorrelation func- 
tion. This supports the idea that the harmonic bath and linear 
coupling approximations are good for this system. 



FMO complex. In Figures 15 and 16, we see that the Har- 
monic approximation gives a roughly temperature indepen- 
dent quantity for all sites. Further, on this scale the spectral 
densities do not differ largely between the different sites. 



2. Harmonic spectral densities for all site of the fmo complex 



Here, we report J(w) the asymmetric component of Giu) 
as described in the text for all temperatures and sites of the 
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Figure 8. Panel a): Two-time correlation function of the energy gap fluctuations of site 1 of the FMO complex, normalized by m'^' at 77K as 
discussed in the text. Panel b): Two-time correlation function for site 1 at 300K. Panel c) Three-time correlation function of the energy gap 

fluctuations of site 1 of the FMO complex, as discussed in the text and multiplied by m'^'/ (^V m(^) j at 77K. Panel d) Three-time correlation 
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Figure 9. Same as for Fig. 8 but for site 2 of the FMO complex. 
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Figure 10. Same as for Fig. 8 but for site 2 of the FMO complex. 




Figure 11. Same as for Fig. 8 but for site 2 of the FMO complex. 




Figure 12. Same as for Fig. 8 but for site 2 of the FMO complex. 




Figure 13. Same as for Fig. 8 but for site 2 of the FMO complex. 
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Figure 14. Same as for Fig. 8 but for site 2 of the FMO complex. 




1200 1500 1800 



Figure 15. Comparison of the asymmetric component of temperature-dependent coupling density Gasym(i^) = J{^) ', for sites 1-4 of the 
Fenna-Matthews-Olson complex, obtained with the Harmonic approximation (As described in the text) at 77K and at 300K. Note that, as 
described in the text, the spectral density can be obtained by dividing J(aj) by 2-k. 
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Figure 16. Comparison of the asymmetric component of temperature-dependent coupling density G'asym(i^) = J{^) ', for sites 5-7 of the 
Fenna-Matthews-Olson complex, obtained with the Harmonic approximation (As described in the text) at 77K and at 300K. Note that, as 
described in the text, the spectral density can be obtained by dividing J{uj) by 2-k. 



